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PURPOSE 

This grant supported research on gas (air) wave bearings for turbomachinery 
applications. A team of researchers from the University of Toledo performed the 
research. The team consists of Dr. Theo Keith (Distinguished University 
Professor), Dr. Florin Dimofte (Senior Research Associate) and Sorin Cioc (PhD 
Research Assistant). 

INTRODUCTION 

The wave bearing has been pioneered and developed by Dr. Dimofte over the 
past several years. This bearing will be the main focus of this research. It is 
believed that the wave bearing offers a number of advantages over the foil 
bearing, which is the bearing that NASA is currently pursuing for 
turbomachinery applications. 

The wave bearing is basically a journal bearing whose film thickness varies 
around the circumference approximately sinusoidally, with usually 3 or 4 waves. 
Being a rigid geometry bearing, it provides precise control of shaft centerlines. 
The wave profile also provides good load capacity and makes the bearing very 
stable. Manufacturing techniques have been devised that should allow the 
production of wave bearings almost as cheaply as conventional full-circular 
bearings. 

RESEARCH ACCOMPLISHMENTS 

The research performed on this grant is described in the following 4 papers: 

1. "Measured Gas Wave Journal Bearing Load Capacity at Steady Loads," (co-authors: 
F. Dimofte, and D.P. Fleming), presented at the 9 th International Symposium on 
Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC-9), Honolulu, 
Hawaii, February 10-14, 2002 


2. "Computation of Pressurized Gas Bearings Using CE/SE Method," (co-authors: S. 
Cioc, F. Dimofte, and D. P. Fleming) STLE Tribology Transactions Vol. 46, No. 1, 
January 2003, pp 128-133. 

3. “Application of the CE/SE Method to Wave Journal Bearings,” (co-authors: S. Cioc, 
and F. Dimofte) STLE Tribology Transactions Vol. 46, No. 2, April 2003, pp. 179-186. 

4. "Calculation of the Flow in High Speed Gas Bearings Including Inertial Effects Using 
the CE/SE Method," (co-authors: S. Cioc, F. Dimofte, and D. P. Fleming) presented at 
the 58th STLE Annual Meeting in New York, New York, April 28 - May 1, 2003 

A copy of each of the first three papers is included below. Because the last paper is under 
review for publication in the STLE Tribology Transactions it is not included in this 
report. 


9 th International Symposium on Transportation Phenomena 
and Dynamics of Rotating Machinery 
Honolulu, Hawaii, February 10-14, 2002 


MEASURED GAS WAVE JOURNAL BEARING LOAD CAPACITY UNDER STEADY LOADS. 


Florin Dimofte, The University of Toledo, currently working at NASA Glenn Research Center, Cleveland, Ohio, 
David P. Fleming, NASA Glenn Research Center, Cleveland, Ohio, 

Theo G. Keith, Jr., The University of Toledo. 


ABSTRACT 

A 35 mm diameter by 28 mm long wave air bearing was 
tested at speeds up to 30,000 rpm. The tests were performed 
on a gas bearing rig at NASA Glenn Research Center in 
Cleveland, Ohio. A maximum load of 196 N was supported 
by the bearing at maximum speed for more than 90 minutes. 
The bearing was stable both dynamically and thermally. Good 
agreement was found with numerical prediction. 

INTRODUCTION 

A circular (plain) journal bearing has better load capacity 
than all other hydrodynamic journal bearings if the 
comparison criteria is the minimum film thickness. However, 
the plain journal bearing is easily de-stabilized, especially if 
the bearing is unloaded. This has been proven both 
theoretically and experimentally [1, 2, 3]. To improve 
stability, the plain journal bearing’s geometry must be 
modified. 

Numerous concepts including lobed, grooved, stepped, 
and tilting-pad bearings were developed. These concepts, 
however, have substantially lower load capacity. The wave 
bearing, a new concept developed at NASA Glenn Research 
Center, offers good stability with a load capacity that is close 
to that of a plain journal bearing [4], The wave journal 
bearing features a wave profile circumscribed on the inner 
bearing diameter. The wave bearing can have 2, 3 or more 
waves, and the wave amplitude can vary. However a three 
wave bearing was found to better satisfy both steady-state and 
dynamic running conditions than other wave bearing 
geometries [5, 6]. A three- wave bearing concept is shown in 
fig. la. A plain journal bearing is shown in fig. Ib for 
comparison. The wave amplitude and bearing clearance are 
exaggerated so that the wave profile can be seen. 

Predictive analysis was developed to quantify the 
performance of wave journal bearings and contrast it to plain 


journal bearings [4]. The predictions revealed a significant 
difference in the pressure distribution between the three-wave and 
the plain journal bearing. The altered pressure distribution 
provides a higher load capacity for the wave bearing when 
compared to the plain circular bearing operating at the same 
eccentricity [4]. In addition, the wave journal bearing offers better 
stability. The analysis was performed assuming that the bearings 
were operating in air, a compressible lubricant. This analysis was 
validated by the measurements performed on a bench test rig that 
was assembled at NASA Glenn Research Center [7], The work 
presented herein was performed to establish the practical maximum 
load capacity of a wave bearing that could also be compared to 
other types of air bearings. 


An air bearing test rig (fig. 2) was built to test journal 
bearings with 30-50 mm diameters and 25-60 mm lengths. A 
commercial air spindle (fig. 2, pos. 4) with a removable add-on 
shaft is used to drive the test bearing. It is mounted vertically to 
eliminate gravitational effects. The spindle can operate at 
rotational speeds to 30,000 rpm with a shaft run-out of less than 40 
micro-inches (1 micron). A cross section of the bearing assembly is 
presented in fig. 3. The mechanical run-out of the add-on shaft 
(fig. 3, pos 4) was measured with a contact electronic dial 
indicator, accurate to 0.1 pm. The run-out of the add-on shaft is 
less than 0. 2 pm at the top (farthest from the spindle) and less 
than 0. 1 pm at the bottom (closest to the spindle). 

The test bearing set up can be seen in fig. 4. The test bearing 
housing (fig. 3, pos. 1 or fig. 4, pos. 1 ) is supported between two 
air thrust plates (fig 4, pos.2 and 3). The thrust plates allow the test 
bearing to move freely in the radial direction. The bottom thrust 
plate (fig. 4, pos. 2) was designed with 3 thrust pads to allow 
bearing alignment to the shaft. Each pad has a separate air supply 
system including a valve (regulator) and pressure gauge. The top 
thrust plate (fig. 4, pos.3) contains a single 360 degree pad. It is 
supported by three threaded rods (fig. 4, pos. 4) which allow easy 
adjustment to the axial clearance of both (bottom and top) thrust 


DESCRIPTION OF TEST RIG 
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plates. 


CONCLUDING REMARKS 


The bearing rotor (fig. 3, pos. 5) is mounted on the 
top of the add-on shaft. The bearing sleeve is assembled 
inside the bearing sleeve assembly (fig. 3, pos. 2). The 
bearing sleeve assembly is fixed inside the housing (fig. 3, 
pos. 1). Rubber O-rings are used between the bearing sleeve 
assembly and the housing to provide additional self- 
alignment capability when load is applied. The test bearing 
was operated in ambient air. 

A radial load can be applied to the test bearing via a 
pneumatic load cylinder (fig.2, pos. 2). A micro-switch (fig. 

4, pos. 5) is used to sense excessive torque on the bearing 
housing; it will shut off the spindle if triggered. A burst 
shield surrounds the experimental bearing for added safety. 

The test rig includes instrumentation to measure shaft 
speed, radial load, rotor-sleeve position, and bearing 
temperature. Shaft speed is measured by a capacitive probe 
located within the spindle. A miniature precision load ceil (fig. 
2, pos. 3 or fig. 4, pos. 6) measures the radial load applied by 
the pneumatic load cylinder. Two light beam displacement 
sensors (fig. 4, pos. 7) allow measurement of rotor-sleeve 
position. The bearing housing temperature is monitored with 
K-type thermocouples (fig. 4, pos. 8). 

TEST RESULTS 

The tested bearing has a 35 mm diameter and is 28 mm 
long. The radial clearance was 18 \im and the ratio of the 
wave amplitude to radial clearance is 0.143. Since the goal of 
this work was to establish a maximum practical load capacity 
of the wave journal bearing, the bearing was tested under 
relatively large loads. Consequently, large eccentricity to 
radial clearance ratios (greater than 0.7) resulted. Speeds of 
5000, 10,000, 20,000, and 30,000 rpm were used. 
Eccentricity, shaft speed, and load were recorded at each load 
step by using a digital camera. A typical photograph 
containing this information can be seen in fig. 5. This figure 
also shows the maximum load applied to the bearing in this 
test, 44 lb (196 N), at a maximum shaft speed of 29,418 rpm 
where the bearing operated with 15.8 |im eccentricity. The 
bearing temperature was also monitored. At maximum speed 
the bearing temperature stabilized at 33°C in a few minutes 
and was then virtually constant throughout the tests, which 
ranged up to 90 minutes. 

The results of this test are plotted in fig. 6 in terms of 
specific load (average bearing pressure) versus eccentricity 
ratio. Analytical predictions are also plotted. Figure 6 shows 
that the wave bearing ran at eccentricity ratios from 0.7 to 

0.88 The bearing run stably at all speeds, and bearing 
temperature stabilized in a short time at reasonable values for 
each run. The wave bearing can carry a specific load of 2 bars 
at 30,000 rpm. 


A 35 mm diameter by 28 mm long wave air bearing was tested 
at speeds up to 30,000 rpm in the hydrodynamic regime. The test 
was performed on the gas bearing rig at NASA Glenn Research 
Center in Cleveland, Ohio. 

A maximum load of 196 N was carried by the bearing at 
maximum speed (30,000 rpm) for more than 90 minutes. This 
corresponds to a specific load of 2 bars. The bearing was stable 
both dynamically and thermally. 


This work was performed at NASA Glenn Research Center in 
Cleveland, Ohio as part of a Director’s Discretionary Fund project. 
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Figure 1. Journal Bearings. 


Figure 2. Gas Wave Bearing Rig: 1 Test Bearing; 2 Pneumatic Load Cylinder; 3 Load Cell; 4 Air Spindle; 5 
Light Beam Displacement Sensor Signal Conditioners. 
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Figure 3. Bearing Assembly Cross Section 
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Figure 5. Records at Max Speed and Max Load: 1 Eccentricity (10|im/div); 2 Shaft Speed (rpm); 3 Load (lb) 
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Figure 6. Specific Loads vs. Eccentricity Ratios: Comparison of Measurement to Prediction 
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TRIBOLOGY-TRANSACTIONS 


Computation of Pressurized Gas Bearings Using 

the CE/SE Method® 


SORIN CIOC, FLORIN DIMOFTE and THEO G. KEITH, JR. 
The University of Toledo 
Department of Mechanical Engineering 
Toledo, Ohio 
and 

DAVID P. FLEMING 
NASA Glenn Research Center 
Cleveland, Ohio 


A numerical scheme that has been successfully used to solve a 
wide variety of compressible flo w problems, including flows with 
large and small discontinuities f entitled the space-time conserva- 
tion element and solution element (CEISE) method, is extended to 

Presented at the 57th Annual Meeting 
in Houston, Texas 
May 19-23, 2002 
Final manuscript approved September 5, 2002 
Review led by R. Gordon Kirk 


compute compressible viscous flows in pressurized thin fluid films. 
This method is applied to calculate the pressure distribution in a 
hybrid gas journal bearing. The formulation of the problem is pre- 
sented, including the modeling of the feeding system. The numer- 
ical results obtained are compared with experimental data. Good 
agreement between the computed results and the test data were 
obtained, and thus, validate the CEISE method to solve such prob- 
lems. 


Nomenclature 

a - radius of the feeding orifice 

a F ,b F ,c F - coefficients in the Taylor series expression of function/ 
a c ,b G> c G = coefficients in the Taylor series expression of function g 
C = radial clearance 

D = diameter of the feeding pocket 

e - eccentricity 

F = fi + gk Vector flux term 

/ = circumferential flux term, see Eq. [5] 

g - axial flux term, see Eq. [5] 

h = film thickness (dimensional) 

h = hIC non-dimensional film thickness 

i,k - unit vectors in circumferential and in axial directions, 

respectively 

L - length of the bearing 

n = Adiabatic coefficient 

n w = number of waves (for wave bearings) 

ns = number of feeding holes 

p - fluid pressure 

P = plp 0 non-dimensional pressure 

p 0 = atmospheric (reference) pressure 

p s = supply pressure 

Q - mass flow rate through one feeding hole 

Q = non-dimensional mass flow rate through one feeding hole - 

see Eq. [17] 

R = journal bearing radius 

t = time 

t 2a ( Q ) 2 non-dimensional time 


u = hp dependent variable in the governing equation 

U - uR velocity in circumferential direction 

U = non-dimensional circumferential velocity 

x = circumferential coordinate 

x = x/(2i tR) non-dimensional circumferential coordinate 

x Q - position of the maximum fluid film thickness, measured in 

the negative direction of axis x 

x w = position of the wave peak (for wave bearings), measured in 

the negative direction of axis x 
z = axial coordinate 

z = z/(2nR) non-dimensional axial coordinate 

y = Isentropic exponent 

A t - time step 

e = weight parameter characterizing one form of artificial 

dissipation, or elC eccentricity ratio 
e w = wave amplitude (for wave bearings) 

ju = fluid viscosity 

p = fluid density 

p = p/pQ non-dimensional fluid density 

p 0 = atmospheric (reference) density 

0 ) = angular velocity of the journal bearing 

Subscripts and superscripts, other than shown above 

( ) 0 = value of the variable at the point O 

( ). = supply pocket index 

( ) x y( ) 2 »( ), = partial derivative with respect to x, z, and t , respectively 
( )" - time step 
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Fig. 1 — Triangular mesh element and its neighbors. 


KEY WORDS 

Bearings; Gas; Hydrostatic Lubrication; Flow Rate 

INTRODUCTION 

Gas lubrication, though preceded by occasional experimental 
work since the mid 19^ centuiy, experienced a strong develop- 
ment in the “Golden Era” of gas lubrication, which started in the 
last years of World War II, and ended in the first part of the 70’s 
(12). In this period much important work was published: (9)-(ll J, 
(2), (6). More recently, significant advancements were reported in 
gas film modeling including rarefaction effects, (13), in numerical 
methods applicable to high-speed bearings, (8), and particularly in 
die treatment of complex geometries, including discontinuities, 

Uh 

The space-time conservation element and solution element 
(CE/SE) method was proposed by (3). Over the past several years 
it has been utilized in a number of fluid flow applications that 
involve shock waves, contact discontinuities, acoustic waves, vor- 
tices and chemical reactions. One of its main features is that it can 
simultaneously capture small and large discontinuities (such as 
sound waves and shock waves) without introducing numerical 
oscillations into the solution, as shown by (5), (4). 

Compressible viscous flow in pressurized thin fluid films, with 
application in hybrid gas bearings, can encounter large pressure 
gradients due to the feeding system or to die large peripheral 
velocities of the bearing. In these conditions, the computational 
methods based on standard finite difference methods or classic 
finite volume methods may prove inadequate having convergence 
problems or inducing oscillations into the solution. Thus, a 
method that is conceptually simple, is second order accurate for 
the entire domain and is able to naturally, deal with large gradients 
and/or discontinuities in the solutions without introducing numer- 
ical oscillations or smearing, is welcome. 


ANALYSIS 

The two-dimensional, transient, Reynolds equation, written for 
a Newtonian compressible fluid in laminar flow is. 



The flow' is considered polytropic, i.e.. 


P 

— = const . , 

p n 


[2] 


where the poly tropic exponent is n = 1 for isothermal flow, n s= 
1.405 for adiabatic flow, or can have other values for general poly- 
tropic flow's. 

A more suitable form of the Reynolds equation for numerical 
formulation is obtained using a new variable u that is the product 
of the non-dimensional film thickness and the pressure, i.e.. 


u~hp 


[ 3 ] 


In terms of w, in non-dimensional variables, the Reynolds 
equation can be written as 


du df dg 
di + dx + dz 


0 


[4] 


where the flux terms / and g are, 



4^—1^ uh *1* 


'fly? 1 — — 

8 “'48 151 

All partial derivatives considered in Eqs. [4] and [5] are carried 
out relative to non-dimensional variables t , x , z . In the following, 
in order to simplify the expressions, the non-dimensional notation 
(upper bar) will be dropped, and all variables wall be implicitly 
considered in non-dimensional form. 

Consider a triangular mesh that covers the (x, z) spatial 
domain. One triangle BCD and its three neighboring elements are 
shown in Fig. 1 . Point A is the centroid of the triangle BCD , while 
points £, F and G are the centroids of the neighboring triangles 
BCH , CD1 and BDJ , respectively. The CE/SE method calculates 
the values of the dependent variables u, u x , u 2 for point A at the 
time step t = t n+ 1 using the corresponding values of the same 
variables for the points £, F and G at the time step t = t \ In order 
to calculate the three unknowns at the new' time step, a system of 
three equations will be derived. 

Consider the quadrilateral ABEC . Simultaneously integrating 
Eq. [4] over the surface of tills quadrilateral and in time, between 
time steps f and t n+ 4 (see Fig. 2), yields 

I f ABEC It- fr dtdu + It- +2 / IaB EC iff + ^) dadt - °* [6] 

Performing the time integration for the first term and trans- 
forming die surface integration into a contour integration for the 
second term (using Green’s theorem) produces 
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Fig. 2 — Conservation volume in the (x, z, /) space. 


I IaBEC ( U " +i - un ) da + ft- +i §A BBC F ■ ndsdt = °> [7] 

where n is the unit vector normal to the contour, oriented outwards 
and 


F = fi + gk [ 8 ] 

is the vector in the (. x , z) plane characterized by the Cartesian unit 
vectors ( i , k). Functions /and g are given by Eq. [5]. Equation [7] 
implies conservation of flux in the three-dimensional space (. x , z, 
t). Functions u, /, g are next written with linear approximations 
using first order Taylor expansions, i.e., 


« - u o + (u*) 0 (s “ x o) + ( u*)o( z ~ z o) + (ttt) 0 (t - <o), [9] 

f-fo + (§£)o( u ~ uo) + (Jf;)o[ u a _ (^x)ol + (5^)0 
\v*z (^z)g] = O'F'U' *F bpU x "h Cp j [10] 


9 ~ 90 + (|^)o( U “ U o) + (^) 0 I u x ~ («*)()] + (5^)0 

fa* - (ti,) 0 l = a GU + + C G - [11] 

Substituting Eq. [9] into Eqs. [10] and [11] yields linear 
expressions for /and g as functions of (x, z, t). In Eqs. [10] and 
[11], coefficients a F , b F , C F , a c , b Q , c G are considered constant 
when integrating Eq. [7], and are known as functions of 
u Qi {u.jc) 0 , (it-) 0 . In Eq. [9], the time derivative can be calculated as 
function of the space derivatives using Eq. [3], i.e., 

( u t ) 0 = ~{fx) 0 — ( Qz)q = — clp(u x )q ~~ O'Gi'U'z) 0 * [12] 

Equations [8]-[i2] are then substituted into Eq. [7], Point 
(X A ',Z A ', t n+ i) is used as the Taylor expansion point for the 
expressions of u 7i+ 5 ,/and g on the contour segments AB and CA , 


while point (%-. Ze> , £ n ) is used as the Taylor expansion point for 
the expressions of //,/and g on the contour segments BE and EC. 
Thus, a first equation with three unknowns, the values «, u z at 
the new half time step is obtained. The coordi- 

nates of points A ' and E’ are selected in a suitable way, as will be 
shown later. The equation is linear and has the general form 

+f>i(u x )^ ! +ci(k 2 )^ 5 + 

diu n B , + ei (u x ) B , + fi{u z ) B , +gi = 0. [13a] 

Two other similar equations are obtained using the same pro- 
cedure for the conservation elements ACFD and ADGB. These 
equations have the general form 

a 2 u n ^ + b 2 {u x ) 1 £* +c 2 (u z )^T 3 + 

d 2 u + e 2 {u x ) B , + f 2 (u z ) B , +g 2 = 0. [13b] 


a 3 u n A ^ +b 3 (u x ) n A t i +c 3 (u z ) n A tC 

d 3 u £, + e 3 (u x )c, + Mu z )q, + g 3 = 0. [13c] 

The linearized system formed by Eqs. [13a], [13b], and [13c] 
can be solved using an iterative method; note that the coefficients 
a b t , c ( ., i = 1,2,3 are functions of the unknowns 
^ alternate approach is to choose the 
Taylor series expansion point A! as the center of the hexagon 
BECFDG. The other three Taylor series expansion points E\ F\ 
G 7 can be chosen arbitrarily; however, in order to maintain con- 
sistency, they are selected as the centers of the corresponding 
hexagons formed around the neighboring triangular elements. 

Adding Eqs. [13a], [13b], and [13c] yields a new equation that 
represents the flux conservation over the hexagon and over one 
half time step. When point A f is the centroid of the hexagon 
BECFDG , this equation has a simpler form given by 

O'suTriU'A' ^ dilCg/ 4 ~ 6i(u x ) pt + 

diUp, + e 2 {u x )p, + /2(u z )y/ + 

d^UQ/ -f- G^{Ux)q, “h fsi^z^G' 9 sum ~ [14] 

where a sum is given by 

O'sum = AabEC + AaCFD + A ADGB = AbBCFDG [15] 

Equation [14] has only one unknown, and can easily be 
solved explicitly. It is also important to note that all the coeffi- 
cients in Eq. [14] depend only on the geometry (coordinates of the 
points) and the values of the dependent variables at the previous 
half time step so that an iterative method is not needed. After cal- 
culating the value the values of the other two dependent 

variables (u x )^ and (u*)^ 5 can be calculated using any two of 
the Eqs. [13a], [13b], and [13c]. This is also called the a scheme. 
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Q = 


\A+(t £> 2 


Q 


U6] 


where the non-dimensional mass flow Q is a fraction C D of the 
ideal mass flow, 


Q = Cd 



for ^ 

for < % < 1 U7] 


When no restricting orifice is present, i.e., an inherent restric- 
tor, Eq. [16] becomes 


Fig. 3 — Feeding orifice geometry. 



Fig. 4 — Comparison between numerical and experimental results for the 
pressure distribution in a gas journal bearing (2 pictures). 


Note that the a scheme introduces “no significant damping” 
(f5J, ( 4 )), so that some form of artificial dissipation is generally 
necessary. The scheme can be simplified and simultaneously sta- 
bilized by calculating the space derivatives in a different way. The 
scheme thus obtained is called the a - e - a - fl scheme. In this 
scheme, the derivatives are calculated as weighted averages 
between the derivatives calculated from the governing equations, 
as shown above (the a scheme), the derivatives calculated using 2- 
D central difference finite difference formulae (weight parameter 
e) and the derivatives using 2-D side finite differencing (weight 
parameter fi). Parameter a is the power index used in the compu- 
tation of the non-linear weighted average that employs 2-D side 
finite differencing. The complete formulation of the way the 
derivatives are calculated in the a - 1 - a - ft scheme can be found 
in (4). It is important that, for a certain value of the weighting 
parameter s (s - 0.5), thta-e-a-/3 scheme eliminates the neces- 
sity of calculating the space derivatives from the governing equa- 
tions, thus the method becomes purely explicit. 

The feed system, formed by a number of ns orifices with the 
general geometry shown in Fig. 3, can be modeled using a gener- 
ally accepted formula, (JO), ( 11 ), that links the mass flow ratio to 
the feed system geometry and the pressures at the ends of the feed 
system considering both effects of the orifice restrictor and the 
inherent restrictor, i.e.. 


Q = xdhyfpJplQ [18] 

The feeding system is introduced into the bearing computation 
through the boundary conditions. Thus, on each feeding pocket 
contour, the mass flow rates calculated from the bearing equations 
and from the feeding system, respectively, must be equal 

(Q0 bearing (Qi) feeding** * [19] 

Equation [19] represents a nonlinear set of ns equations, where 
. the pressures in the bearing at the feeding holes are the unknowns 
p r i = 1,2 ,...,ns. This system is solved iteratively using Newton’s 
method. 

RESULTS AND DISCUSSION 

A computer code has been developed based on the described 
method, using the a-e-a-t 3 scheme with e = 0.5. This code has 
been tested for some simple cases where experimental data were 
available, considering isothermal flow (n = 1). 

The first case considered was the circular gas bearing without 
any feeding system. The results are shown in Fig. 4, where the 
same bearing with the aspect ratio LI(2R) = 2.0 is examined for 
two bearing numbers A = (6p «/T)/(p 0 C )and for two eccentrici- 
ties. There are no significant differences between the experimen- 
tal and calculated pressure distributions, at least in the middle 
plane where the experimental results are available. The relative 
differences between the experimental and calculated non-dimen- 
sional loads t, = Load/(j? 0 2LR) are 3.5% and 4.0%; both represent 
improvements over the theoretical results shown by (6). 1 

The second case considered is the wave journal bearing (7) 
without any feed system. The fluid film thickness in an aligned 
wave bearing is given as, 

h= 14- ecos2'n{x + a’o) + e w cos[27rn w (x + a;**,)] [20] 

In the last term on the right hand side of Eq. [19], e is the 
wave amplitude. n w is the number of waves and x w is the wave ori- 
gin relative to the origin of the x axis. For the present case, the 
bearing diameter and bearing length are 50 mm, the radial clear- 
ance is 0.02 mm, n w = 3, = 0.3, and the relative minimum film 

thickness is h m j n = 0.3 for all cases. The bearing number is 
3.566. 

The results are compared with a finite difference based code 
built by (7). Different bearing positions Cx lv ) have been tested and 
the results (calculated load magnitude and position) are presented 
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Table 1 — Comparison between Calculated Data with CE/SE Method and with Finite Difference 

Method 

Eccentricity 

a w 

(deg) 

Load CE/SE 
(N) 

Load FD 
(N) 

Load Angle 
CE/SE (deg) 

Load Angle 
FD (deg) 

0.502 

40 

177 

177 

40.49 

40.18 

0.460 

32 

170 

170 

39.09 

38.8 

0.433 

24 

166 

165 

36.37 

36.11 

0.416 

16 

163 

162 

32.94 

32.69 

0.404 

8 

158 

158 

29.30 

29.06 

0.400 

0 

154 

154 

25.61 

25.38 

0.404 

-8 

151 

151 

22.08 

21.86 

0.416 

-16 

149 

148 

18.95 

18.75 


Table 2 — Pressurized Gas Bearing Geometry and Working 
Conditions 

Bearing length (mm) 

117.5 

Bearing diameter (mm) 

60.4 

Supply planes 

2 

Supply plane position (mm) 

12.7 

Holes/supply plane 

14 

Orifice diameter (mm) 

0.16 

Pocket diameter (mm) 

0.9 

Supply pressure (Pa) 

5.514x10 s 

Injection angle (deg) 

90 


Position of Jet 



in Table 1 . The results show very good agreement between the two 
methods. 

The third case considered is a pressurized gas bearing without 
rotation and with zero eccentricity. Details of the bearing geome- 
try and working conditions are presented in Table 2. Figures 5(a) 
and 5(b) show the calculated and experimental pressure distribu- 
tions for two values of the radial clearance, corresponding to the 
subsonic flow regime (obtained when C = 12.7 pm), and the 
choked flow regime (obtained when C = 31.75 pm) in the feed 
system, in two longitudinal planes situated at the jet position (Fig. 
5(a)) and at half distance between jets (Fig. 5(b)); pressure peaks 
are visible at the jet positions; also nearly constant pressure is 
obtained between the supply planes. The predicted pressure peaks 
at the position of the jet are higher than the experimental values 
(this difference is more visible for the subsonic inlet flow, howev- 
er it is present in both cases). This is due to the difficulty of meas- 
uring the local pressure at the feeding orifice position. The value 
for the correction factor C D is 0.8 for the subsonic inlet flow, and 
0.85 for the choked flow. 

Finally, the code was applied to different configurations, where 
experimental or computed data were not available. 

Figure 6 shows the pressure distribution (sections at the jet and 
at the mid-jet positions) obtained for five supply planes for the 
same bearing as in the previous case C = 12.7 pm. The supply 
planes are equally distanced with respect to each other and with 
the bearing ends. It is seen that, although the external pressure is 
the same for all supply holes, the pressure that develops in the 
pockets (inside the bearing) is not the same for all supply planes. 


Fig. 5(a)—Comparison between the calculated and experimental pres- 
sure distributions at jet position. 


Mid Jet Position 



Fig. 5(b)— Comparison between the calculated and experimental pres- 
sure distributions at mid-jet position. 


The pressure distribution between two consecutive feeding planes 
has an almost linear form. Furthermore, at the central supply 
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Fig. 6 — Pressure distribution for a circular bearing without eccentricity, 
with five supply planes. 

plane, the peak pressures are not as visible as the peak pressures 
at the other supply planes; this suggests that the central supply 
plane does not have an important contribution for the general 
pressure distribution inside the bearing. 

Figure 7 shows the pressure distribution in half bearing for the 
same bearing geometry as in the previous two cases, but in a run- 
ning condition (40,000 RPM , A = 25.76), with a relative eccen- 
tricity £ = 0.5. The value C D = 0.8 was used. The pressure distri- 
bution is very complex. It can be seen that for many feeding holes 
the flow is inverted, i.e., is directed from the bearing towards the 
feeding system, because the pressure inside the bearing is higher 
than the supply pressure. 

CONCLUSIONS 

The CE/SE computational method has been extended for the 
first time to calculate compressible viscous flow in pressurized 
thin fluid films with restrictors in series. Formulation of the 
method was presented along with the numerical results obtained 
for both non-press urized and pressurized bearings. The results 
were compared with existing experimental and computed data. 
The results demonstrate the ability of the method to accurately 
predict the pressure distributions in such flows. 

While most numerical methods handle space and time differ- 
encing terms in the governing equations separately, this relatively 
new scheme treats them in a unified way. This results in very good 
accuracy even though the unknowns are considered locally linear 
and even when a relatively course grid is used. Also, the entire 
structure of the method, as well as the developed code, are rela- 
tively simple. Unlike most numerical schemes, no special treat- 
ment is necessary for discontinuities, providing that the governing 
equations are written in strong conservative form. However, the 
method is limited to time dependent equations. Steady state 
results, such as those discussed in this work, can be obtained only 
after the stabilization of die unsteady solution starting from initial 
conditions. 



Fig. 7 — Pressure distribution for a circular bearing with eccentricity, 
with five supply planes. 


Future work will focus on predicting dynamic characteristics 
of pressurized gas bearings. Improving the feeding system with 
restrictors in series modeling will also be considered, as well as 
improving the computation time. Cases including discontinuities 
will also be considered. 
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INTRODUCTION 

The space-time conservation element and solution element 
(CE/SE) method was proposed for the first time by Chang and To 
(]). Over the past several years it has been utilized in a number of 
fluid flow applications that involve shock waves, contact discon- 
tinuities, acoustic waves, vortices and chemical reactions. Being 
capable of simultaneously capturing small and large discontinu- 
ities (such as sound waves and shock waves) without introducing 
numerical oscillations in the solution (2)-(4) this new method is an 
excellent candidate to be applied to the flow in cavitated wave 
bearings. 

Historically, the most common practice to account for the 
effects of cavitation in fluid film bearings was through the appli- 
cation of the Giimbel (or half-Sommerfeld) boundary conditions. 
This was accomplished by setting negative pressures to zero (rel- 


Nomenclature 

a - coefficient in the flux term /(see Eq. [4]), m/s 

b = coefficient in the flux terms / and g (see Eq. [4]), m 2 /s 

c = coefficient m the flux term g (see Eq. |4}), m/s 

C = radial clearance, m 

d m = degree of misalignment, non-dimensional 

f = Flux term in x direction (see Eq. [3]), m 2 /s 

F - fi 4- gk vector flux term 

8 = flux term in z direction (see Eq. [3]), m 2 /s 

g { , = switch function (see Eq. [5]), non-dimensional 

h = film thickness, m 

ifk = unit vectors in circumferential and in axial directions, 

respectively. 

L = bearing length, m 

n = unit vector normal on the contour line, oriented outwards 

P ~ fluid pressure, Pa 

R = bearing radius, m 

s = linear (contour) coordinate, m 


t = time, s 

u - /?6, m 

U = velocity in circumferential direction (a )R), m/s 

x - circumferential coordinate, m 

i = axial coordinate, m 

a = angle between the centerline and the direction of the mis- 

alignment at the bearing center 
p = bulk modulus, Pa 

A t = time step, s 

)t = fluid viscosity, Pa*s 

0 = p/p r non-dimensional density in full film: fractional film 

content in cavitated region 

cu = angular velocity of the journal bearing, s' 3 

Subscripts and superscripts 
( ) 0 = value of the variable in the point O 

( ) t = value at the cavitation border 

( ) x , ( h 

( ) t = partial derivative with respect to .v, z, and /, respectively. 
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ative to the cavitation pressure). Although the load carrying pre- 
dictions by this approach were reasonably accurate, the results 
violated the mass conservation principle. Consequently, several 
other procedures have been proposed. Jakobsson and Floberg (5) 
and later Olsson (6) introduced a self-consistent, mass conserva- 
tive, set of boundary conditions for cavitation to be applied to 
Reynolds equation. This procedure is valid for moderately to 
heavily loaded bearings and is generally called JFO theory. This 
methodology is commonly incorporated into modern computa- 
tional algorithms for bearings, and is also implicitly included in 
the present method. 

Previous computational methods used for this problem are 
known to have certain difficulties. Elrod’s algorithm (7) was the 
First that correctly applied JFO theory, however the algorithm 
necessitated, as the author pointed out, “considerable experimen- 
tation” to develop. Moreover, it has only first order accuracy in 
the cavitated region, while in the full film region the algorithm is 
accurate to second order; an oscillation in the cavitation front is 
often found to occur. The method proposed in (8) is based on con- 
cepts used in transonic flow (9). This method uses a number of 
features from the Elrod algorithm, but does not “rely on experi- 
mentation” to develop the solver. It should be noted that the 
Vijayaraghavan and Keith method has the same accuracy as the 
Elrod’s algorithm in cavitated regions and in the full film region. 
In addition, like Elrod’s method, the solver loses accuracy at the 
cavitation boundaries. 

In this context, a method such as the CE/SE scheme, which is 
conceptually simple and has the capacity to accurately predict the 
fluid film flow including the boundaries of the cavitated region(s), 
without numerical oscillations or smearing, represents an 
improvement over these earlier methods. This is especially impor- 
tant since, as shown in (10) for the one-dimensional case, the the- 
oretical formulation can lead to flow discontinuities at the refor- 
mation front. The CE/SE method applied to the Reynolds equation 
has some noticeable advantages: it is second order accurate over 
the entire domain, it computes in a unified way the pressure 
induced flow for all the regions and, because it solves a set of inte- 
gral equations derived directly from the physical conservation 
laws, the scheme is able to capture naturally the flow discontinu- 
ities (cavitation boundaries). Thus, the CE/SE method can poten- 
tially provide better results for the cavitation and reformation 
front positions, as well as for the distributions of the state vari- 
ables in their vicinity. This paper presents the general develop- 
ment of the CE/SE method applied to two-dimensional cavitated 
flows in journal bearings, showing the strengths and some of the 
limitations of this emerging modem computational scheme. The 
paper also presents some results obtained using this numerical 
scheme to illustrate its characteristics. 


ANALYSIS 


The theoretical model, based on Reynolds equation in Elrod’s 
formulation, is the two-dimensional extension of the formulation 
presented in (10). In strong conservation form, required for a 
weak solution (with discontinuities), the governing equation is, 


du Of dg 

1 q - 0 

dt dx dz 


[1] 


The main unknown u is the piouuci of ihe fiim thickness, h, 
and the non-dimensional density, 8, 


u = /i0 


[ 2 ] 


Note that the non-dimensional density 6 has also the meaning 
of fractional film content in the cavitated region. The flux terms 
are given by the following expressions, 

du 


f = au — b 


dx 


g = cu — b 


du 

dz 


where coefficients a, b, and c are 

_U j3h Oh 
° - 2 + Qc 12 j_i dx 


9c 


(3h 2 ■ 
12m 


(3 h dh 

C ~ 9c \2 /x z 


[3] 


[4] 


9c 


[5] 


The value of the switch function g c is determined by the value 
of 8, 

1 — full film region , 9 > 1 
0 — cavitated region 6 < 1 

The governing equation is spatially elliptic in the full film 
region and hyperbolic in the cavitated region. 

In order to numerically solve the problem using the CE/SE 
method, the main unknown u and the flux terms /and g are writ- 
ten in linear form using a Taylor series. Considering a generic 
point of expansion O, they are written as 


u = uq + {u x )o{x - s 0 ) + (u z ) 0 (z - z 0 ) 4- (u t )o[t - to) [6] 

/ — fo + (fx)o{x — ^o) + (fz)o{ z “ Zq) + {ft) o(£ - <o) 

(J — 9o + (f/a:)o(-T ” ^’o) + (Oz)o{ z “ z o) + {(0t)o{t ~ Ui) V] 

Because the time derivative u t can be written from the govern- 
ing equation in terms of space derivatives, Eq. [6] can also be 
expressed as 


u = Uo + ('U;,)o('C - So) + (u z )o{z - Zo) - ( fa ; + fj z )o{t - t 0 ) [8] 

Solution of the governing equation, Eq. [1], is performed on a 
triangular mesh in the (x, z) plane. One triangle BCD and its three 
neighbor elements are shown in Fig. 1. Point A is the centroid of 
the triangle BCD , while points E , F and G are the centroids of the 
neighboring triangles BCH , CDl and BDJ. The CE/SE method 
calculates the values of the variable w, for a “representative point” 
A’ of the triangle BCD, at the new time step / « r n+l using the cor- 
responding values of the same variables for the points E\ F y and 
G ' at the old time step t = t n . The location of the “representative 
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Fig. 1— Triangular mesh element and its neighbors. 


points” A\E\F\ and G’ will be explained later. In order to cal- 
culate die unknown at the new time step, one equation will be 
derived. This equation can be obtained considering the hexagon 
BECFDG. Simultaneously integrating the governing equation, 
Eq. [1], over the surface of this hexagon and in time, between time 
steps f and r” +3 , yields 


[ f f —dtdxdz + f f f + ^)dxdzdt = ( 

J BE JCFD Ji n oi Jbe Jcfd dz oz 


[9] 


Equation [9] implies flux conservation in the three-dimension- 
al space (x, z, r), which represents the central feature of the 
method, also justifying its name. Performing the time integration 
for the first term and transforming the surface integration into a 
contour integration for the second term (using Green’s theorem) 
yields, 


JbeJcfd^’ 1 * 1 ~ u ’ l ) dxdz + fF §becfd ? ' ndsdt = 0 I'O] 


where ri is the outward directed unit vector, normal to the contour, 
and 


F = fi + gk [11] 

is a vector in the (x, z) plane characterized by the Cartesian unit 
vectors ( i . k). 

Because functions w, /, and g are substituted with linear 
approximations, given by Eqs. [7] and [8], the time integral in Eq. 
[10] can be written as the value of the integrand at the mid-inter- 
val multiplied with the length of the time step, i.e., 

1 BE JcFd( U ’ 1+1 - + WfBBCFV P • «^) n+,/2 = 0 H2] 

The next important feature of the method consists in selecting 
the “representative point” for each triangular element as the cen- 
troid of the hexagon constructed around that triangular element. In 
this case, point A’ is the centroid of hexagon BECFD , while points 
E\F\ and G’ are the centroids of the hexagons corresponding to 
the neighbor elements. 


The integrals that appear in Eq. [12] are calculated using the 
linear approximations, given by Eqs. [7], and [8]. The points of 
expansion are as follows: 

- Point A’ for the first term in Eq. [12], which is approximated 
as, 

/ / u n+1 da = u'^Abecfd 113] 

J BE JCFD 

- Point E* for all the integrals that include the quadrilateral 
ABEC, and similarly, point F’ for the quadrilateral ACFD , and 
point G ’ for the quadrilateral ADGB. The surface integral 
Jbe Jcfd un d& is written as the sum of the surface integrals 
over the three mentioned quadrilaterals, each of them with its own 
expansion point, thus, 

Ibe Jcfd undo “ Sab Sec v " da + Sac Sfd + Sad Jcb uUd<T f * ^] 

For example, the first integral in the right hand side of Eq. [14] 
is written as, 


IaeJec u " da - Iab Iec I"e' + (“*)*(* - *£') + - ZB')]dxdz[ 15] 


Since all the terms are evaluated at the old time step n, the right 
hand side can be determined numerically. Similarly, the contour 
integral § BECFD F • fids is written as the sum of six segment 
integrals. 



J BECFD 


fids = 


(/ F-nds+ I F ■ nds)E ,J r 
Jbe Jec 


(/ F- fids A- / F-nd$)p'+ 

Jcf Jfd 


([ F* fids A- [ F-nds)c' [16] 

Jdg Jgb 

The right hand side of Eq. [ 16] was written in segregated form 
to indicate that each group of terms uses one expansion point. For 
example, point £’ is used as expansion point for the integrals over 
the segments BE and EC, etc. All the line integrals are calculated 
using the values at the old time step, which are readily available. 
For example, 

( / F ■ ?ida ) n+1/2 = 

Jbe 


/ I F%, + (F?) EI <x-x E ')+ 

Jbe 

_ _ At 

(F:)e’ + (■ F?)e ' (2 - ZB') + ( F?)b ' -y] • nda [1 7] 

At the end of this process, one explicit equation is obtained. It 
contains only one unknown, viz., the value of u ,7+] at the new time 
step. 

After calculating the field of the unknown u at the new time 
step, the space derivatives u x and u, must also be calculated. This 
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Fig. 2— Uniform triangular mesh. 


task is accomplished by using an averaging combination of cen- 
tral and side differencing approximations. This procedure insures 
a minimum, but necessary, amount of numerical dissipation. For 
more information about the calculation of the space derivatives, 
the reader should consult (11) and (12). 

An important feature is that, even though the space derivatives 
are not calculated from the conservation condition (the integrated 
form of Reynolds equation), this is still satisfied over the surface 
of each hexagonal element and in time (conservation volume). 

The flux terms derivatives f x , / r , g x> g 7 can be calculated fol- 
lowing an identical procedure. Another approach, which was also 
used, considers, starting from Eq. [3], and neglecting higher order 
terms, so that 

(/x)c = (u*)o, (Jz) o - oo(*)o.(/t)o - Oo(u<)o = -ag(u s )o - aoCb(uz)o 
( 9 x)o - co(wx)o, (52)0 = co(n.)o,te)o = Co(u t )o = -aoCo(Ur)o - Co(u = )o [ 18 ] 

These approximations have the advantage of simplifying the 
calculations, without significantly reducing the accuracy of the 
method. 

An interesting feature of the method consists in the fact that the 
value of the unknown at the new time step for any element is cal- 
culated using exclusively the values at the old time step for the 
neighboring elements. Figure 2 shows a portion of simple uniform 
triangular mesh. Assuming that the values of u at the time step n 
are known for all elements marked with black dots, at the time 
n+ 1, the values of u for all elements marked with hollow circles 
can be calculated. At the next time step, n+ 2, the values for the 
black dot elements can be calculated. It is thus evident that the 
method actually uses two staggered grids, and at each time step it 
switches from one grid to the other. In practice, both grids can be 
used at every time step. This procedure doubles the volume of cal- 
culations, but it has the advantages of producing better resolution 
of the results, as well as permitting the extension of the applica- 
bility of the method for general unstructured grids. Another 
advantage is related to the calculation of the space derivatives, as 
shown in (12). 

As with most explicit methods applied to hyperbolic partial 
differential equations, the stability of the method is subject to the 


CFL condition ( 13), which states that the ratio between the time 
step and the space step must be small enough (more exact, the 
maximum Courant number, which is the product of this ratio and 
the local velocity, must be smaller than 1). Also, when the time 
step (Courant number) decreases, the accuracy of the method 
decreases, because the dissipation increases. For the present prob- 
lem, the Courant number is multi-dimensional and variable in the 
field. In addition, the governing equation is second order, necessi- 
tating the calculation of an equivalent Courant number. Because 
of these features, a different approach is used. In this methodolo- 
gy, the time step is determined using a heuristic procedure, so that 
it simultaneously satisfies the following conditions: (i) when the 
time step is increased, by a factor of two for instance, the scheme 
becomes unstable, and (ii) when the time step is decreased, by a 
factor of two for example, the results do not change significantly 
(the change is due only to a small increase of the dissipation - 
smearing of the solution). This approach is directly derived from 
the general characteristics of explicit numerical methods, and it 
was found to be very effective. An example of a detailed approach 
for the convergence and error bound analysis of this method 
applied at a more simple equation can be found in (76). The use of 
a variable time step did not bring any significant benefit. 

RESULTS AND DISCUSSION 

In order to determine the performance of the method, numeri- 
cal solutions in several geometries were performed. The results 
were compared with the results obtained using other numerical 
methods. Only the asymptotic steady state solution, obtained by 
time integration until the state parameters stabilize, was consid- 
ered. The grids used consisted of 90-180 intervals in the circum- 
ferential direction (a courser grid can be used for simple geome- 
tries), 10-20 intervals in the axial direction for a half bearing, 20- 
40 intervals for a full bearing (misaligned cases). The grids used 
are uniform in both directions, and the results show no significant 
change for the grid intervals specified. The non-dimensional time 
step t = cor has the order of 10‘ 5 . 30,000-120,000 time steps were 
usually needed to obtain the steady solution (residual of variable 
u smaller than 10" 4 - 10' 5 ). The computational time required to 
obtain the steady-state solution starting from an initial uniform 
pressure distribution ranged from 5 to 25 minutes on a 1 GHz PC. 

Circular Journal Bearing 

A standard journal bearing with one inlet groove was first con- 
sidered, The geometry and the fluid characteristics are presented 
in Table 1. 

Figures 3 and 4 show a comparison between the non-dimen- 
sional gage pressure, p = and the fractional fluid film 

content, 0, distributions obtained with the CE/SE method and with 
the type difference method (9), (14). Two transverse sections 
through the bearing at z = 0 and z = L/4 are presented for each 
method. The bearing z coordinate ranges between the values -L/2 
and -L/2, so that z = 0 is the symmetry plane. The differences 
between the two methods are seen to be relatively small; the max- 
imum pressure predicted using type differencing is less than 2% 
lower than the maximum pressure calculated with the present 
method, while the total load is 6.9% lower. The difference 
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Table 1 — Physical Conditions for Circular Journal Bearing 

Parameter 

Value 

Units 

Length 

62.8 x 10' 3 

m 

Diameter 

62.8 x 10' 3 

m 

Clearance 

40.0 x 10' 6 

m 

Relative eccentricity 

0.8 

- 

Groove position 

0 

deg 

Supply pressure (gage) 

1.355 x 10 6 

Pa 

Angular velocity, oj 

628.32 

s' 

Viscosity 

0.0035 

Pa s 

Bulk Modulus, p 

5.421 x 10 7 

Pa 

Cavitation pressure (gage) 

0 

Pa 



Fig. 3 — Pressure distribution for a circular journal bearing with one inlet 
groove. 


between the attitude angles calculated with the two methods is 
0.9*. 

Figure 4 shows that for the cavitated region, 6 < 1 , the type dif- 
ferencing method produces a sharper discontinuity in the fraction- 
al film content compared to the CE/SE method. However, in the 
type differencing code, the discontinuity in the fractional film 
content distribution is forced to occur using different numerical 
formulations for tine two different regions - full film, and cavita- 
tion, respectively. The computational domain is accordingly split 
into two parts, connected through boundary conditions. On the 
other hand, the computational domain used by CE/SE code is con- 
tinuous (periodic boundary conditions are used at one circumfer- 
ential position), so that the fractional film content discontinuity' 
appears in the field. Also, no special treatment is used for any of 
the two regions. This approach was possible only because the 
CE/SE method is able to cope with large discontinuities without 
introducing any significant amount of numerical smearing and/or 
oscillations. The code developed using the CE/SE method is thus 
more general and potentially more robust. The visible presence of 
numerical dissipation is due to the small time step, which is nec- 
essary to obtain the CFL stability condition in the full film 
regions. This dissipation affects only the value of 0 in the cavitat- 
ed region in the vicinity of a discontinuity, which is unimportant 



Fig. 4— Fractional film content distribution for a circular journal bearing 
with one inlet. 
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Fig. 5— Pressure distribution in a wave bearing. 


relative to the pressure distribution, since p = p c = const, in these 
regions. 

Wave Journal Bearing 

The geometry of a wave bearing is more complex compared 
with the geometry' of a standard journal bearing, Dimofte (15). An 
example of the non-dimensional film thickness h — hjC distri- 
bution for a three- wave bearing as a function of the circumferen- 
tial coordinate is shown in Fig. 7(b), as discussed later. The phys- 
ical conditions are presented in Table 2. 

Figure 5 shows the non-dimensional gage pressure 
p = ^(^) 2 distribution obtained using the present method 
compared with the pressure calculated by Dimofte, using a finite 
difference (FD) method to solve the steady form of Reynolds 
equation with Giimbel (or half-Sommerfeld) boundary conditions. 
The results show similar variations, however the peak pressures 
differ. The total load predicted by the present method is 2039 N, 
13.7% higher than the load predicted using FD, which indicates a 
lower numerical dissipation for the CE/SE method. The difference 
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Table 2 — Physical Conditions for Wave Bearing 

Parameter 

Value 

Units 

Length 

26.0 x 10' 3 

m 

Diameter 

45.0 x 10' 3 

m 

Clearance 

15.0 x IO' 6 

m 

Number of supply pockets 

3 

- 

Pocket positions 

86, 206, 326 

deg 

Pockets width 

4.0 


Supply pressure (gage) 

5.458 x 10 5 

Pa 

Angular velocity, u) 

413.85 

-l 

s 

Viscosity 

0.00325 

Pa s 

Density 

902.0 


Bulk Modulus, p 

1.2105 x 10 s 
1.0290 x 10 7 

Pa 


Table 3 — Physical Conditions for Misaligned/Aligned Wave 
Bearing 

Parameter 

Value 

Units 

Length 

45.0 x 10' 3 

m 

Diameter 

45.0 x 10' 3 

m 

Clearance 

15.0 x 10' 6 

m 

Angle of misalignment 

90 

deg 

Degree of misalignment 

0.5,0 

- 

Number of supply pockets 

1 

- 

Pocket position 

80 

desL 

Pocket width 

4.0 x 10' 3 

m 

SuddIv pressure (gage) 

5.458 x 10 5 

Pa 

Angular velocity, co 

413.85 

s l 

Viscosity 

0.00325 

Pa s 

Bulk Modulus, ($ 

1.2105 x 10 8 

Pa 


between the predicted load directions using the two codes is less 
than 1.6 . These results have been obtained using a bulk modulus 
value (3 = 1 .2105 • 10 8 Pa. The pressure distribution using a small- 
er value, P = 1.029 • iO 7 Pa, is presented in Fig. 6. In this case, 
even though both peak pressures predicted by the present method 
and FD are almost the same, the total load compared with the pre- 
vious case does not change significantly. The load directions how- 
ever have a larger difference (5.5 ). In brief, the film compress- 
ibility effects reduce the pressure peaks in the bearing and also 
produce a change in the pressure distribution phase. Another 
effect, this time related to the numerical scheme, resides in the fact 
that when the value of the bulk modulus increases, the diffusion 
velocity increases. This leads to a smaller maximum time step for 
the scheme to be stable. The method more readily “handles” com- 
pressible fluids. Note that, although the bulk modulus value for 
the oil used by Dimofte is not known, the first case (p = 1.2105 • 
1 0 8 Pa) probably is a more realistic value. Also, the code based on 
the present method uses a film thickness distribution calculated by 
Dimofte considering the elastic deformation of the bushing; this 
can also be a source of errors. 

Misaligned Wave Journal Bearing 

The misalignment of journal bearings can be defined using two 
parameters: the angle a between the centerline and the direction 
of the misalignment at the bearing center and the degree of mis- 


Pr«*ure {SstfltxitlOfl 



Fig. 6— Pressure distribution in a wave bearing. 

alignment d m that represents the proportion of the actual mis- 
alignment to the maximum possible. The maximum misalignment 
is restricted by the condition that at the bearing ends (in the axial 
direction) the film thickness reaches the value of zero. The physi- 
cal conditions of the bearing are presented in Table 3. Because of 
the misalignment, the fluid film thickness is a function of both the 
circumferential and axial coordinates, so that for each value of the 
circumferential coordinate there is a domain of thickness varia- 
tion, represented in dark color in Fig. 7(a). Both the aligned and 
misaligned bearing have the same film thickness distribution at 
the center plane (which is also the symmetry plane for the aligned 
case). 

Three sections through the pressure distribution, at the mid- 
plane z = 0, and at mid-distance between this plane and the bear- 
ing ends, z = ± L/4, are presented in Fig. 8(a). The maximum cal- 
culated pressure is 38.8 MPa (not seen in the figure). The same 
bearing with identical physical conditions but without misalign- 
ment has different values for the pressure distribution, as seen in 
Fig. 8(b) for the same sections (the pressure distribution is sym- 
metric relative to the central plane z = 0. The maximum calculat- 
ed pressure in this case is 23.2 MPa, at the symmetry plane. The 
cavitation and the full film regions for the misaligned and the 
aligned 3-wave bearing are presented respectively in Figs. 9(a), 
9(b). The load is found to be 10,523 N for the aligned bearing and 
22,115 N for the misaligned bearing, which shows that the mis- 
alignment may have a positive impact on bearing performance. 
This effect is due to the fact that, with misalignment, a region with 
smaller film thickness is present, compared with the same aligned 
bearing. This determines the development of higher pressures. At 
the same time, the region that has a larger film thickness develops 
cavitation, so that the pressure does not decrease the same amount 
as it increases in the smaller films thickness region. 

CONCLUSIONS 

The CE/SE method was applied for the first time to investigate 
two-dimensional flow in cavitated wave bearings. The theoretical 
formulation of the solution method was presented along with 
numerical results. The results were compared with the results 
obtained using other numerical algorithms. Using Elrod’s formu- 
lation, discontinuities can appear at the reformation fronts, and the 
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Fig. 7(a)— Fluid film thickness domain in a misaligned 3-wave bearing. 


Fig. 8{a)— Pressure distribution in a misaligned 3-wave bearing at dif- 
ferent axial sections. 



Fig. 7(b) — Fluid film thickness in an aligned 3-wave bearing. 

results presented show that the CE/SE method is able to capture 
these discontinuities without any special treatment, and thus it is 
potentially more accurate than “standard” methods. It is also 
important to recognize that the CE/SE method is conceptually 
simple and entirely explicit, which makes it also computationally 
efficient. However, the computational time is higher Ilian for the 
steady solvers, which can give the results in seconds on modern 
personal computers, but it is comparable with other transient 
solvers (this statement is based on the comparison with the tran- 
sient type difference code, which also uses the Elrod's formula- 
tion). The conclusion is that the method is a strong candidate to 
solve applications that require more precise results, such as accu- 
rate, robust computation of the cavitation boundaries, as well as to 
solve transient problems. 
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Fig. 9(a) — Cavitation map in a misaligned 3-wave bearing (full film 
region shown in dark color). 
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